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Abstract 

We develop a coarse grained (CG) approach for efficiently simulating calcium dynamics in the 
endoplasmic reticulum membrane based on a fine stochastic lattice gas model. By grouping neigh- 
boring microscopic sites together into CG cells and deriving CG reaction rates using local mean 
field approximation, we perform CG kinetic Monte Carlo (kMC) simulations and find the results 
of CG-kMC simulations are in excellent agreement with that of the microscopic ones. Strikingly, 
there is an appropriate range of coarse proportion m, corresponding to the minimal deviation of 
the phase transition point compared to the microscopic one. For fixed m, the critical point in- 
creases monotonously as the system size increases, especially, there exists scaling law between the 
deviations of the phase transition point and the system size. Moreover, the CG approach pro- 
vides significantly faster Monte Carlo simulations which are easy to implement and are directly 
related to the microscopies, so that one can study the system size effects at the cost of reasonable 
computational time. 
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I. INTRODUCTION 



As a second messenger in living cells, calcium ions (Ca^"*") play a vital role in providing 
the intracellular signaling. Many important cellular processes and biological function, such 
as muscle contraction and synaptic transmission, are regulated by Ca?^ signals 
release is an inherently multi-scale problem, for instance, in cardiac myocytes, the majority 
of the control of calcium-induced-calcium- release (CICR) js, 6| happens in the microdomain 
of the so-called diadic cleft, this microdomain is between the L-type voltage-gated Ca?^ 
channels and the ryanodine receptors. The ryanodine receptors 'sense' local [Ca^"*"] in the 
diadic cleft positioned between the t-tubules and the sarcoplasmic reticulum. The length 

scale of aforementioned occurrences is on the order of nanometers and relevant time scales 

n 



range from micro- to milliseconds 
diadic spaces which act independently 



iowever, each cell contains approximately 10,000 
sl. Typically, one is interested in Ca?^ currents 
at the whole cell level and higher. This is a multi-scale problem. Given the limits of 
computational power, hardly can we model an entire cellular cytoplasm by incorporating 
detailed structural information. 

Some multi-scale models of CICR have been developed that successfully reproduced ex- 
perimental observations, as well as save computation largely [t-^. However, these models 
are based upon deterministic coupled ordinary differential equations derived from biophys- 
ical mechanisms , and lack accurate description of microscopic dynamics of calcium ion 
channels. In fact, fluctuations are always exist in ion channels and play a crucial role in 



Co?^ release mechanism 
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12| . Recently, Vlachos and coworkers proposed a multiscale 



approach for coarse graining stochastic processes and associated Monte Carlo (MC) simula- 



tions in surface reaction systems 



ISMlq. The method is efficient in describing much larger 



length scales than conventional MC simulations while still incorporating microscopic details, 
and resulting in significant computational savings. An overview of the method is given in 

In the present work, the multiscale approach was applied to a relatively simple stochastic 
reaction-diffusion lattice model for calcium dynamics in the endoplasmic reticulum (ER) 



membrane, proposed by Guisoni [17|, |18[. We coarse grain the model and processes, and 
derive the coarse-grained (CG) surface diffusion transition probability rates. By numerical 
simulations, it is found that the results of CG kinetic MC (kMC) simulations are in excellent 
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agreement with that of the microscopic ones corresponding to the optimal coarse proportion. 
Secondly, we study the system size effects by fixing the coarse proportion, and find the phase 
transition point increases monotonously as the system size increases. Especially, there exists 
a scaling law between the deviations of the phase transition point and the system size. 
Finally, we investigate CPU time and find the approach provides significantly faster MC 
simulations which are easy to implement and are directly related to the microscopic one. 



II. COARSE-GRAINING THE LATTICE MODEL 



Microscopic Model — We consider a two-dimensional square lattice with two interpenetrat- 



ing sublattices A and B [17|, ll8|| in ER membrane, as shown in Fig. [T] Calcium channels are 
located only on the sites of the sublattice B and calcium ions occupy not only the sites of 
the sublattice A but also the sites of the sublattice B. A site i of the sublattice A can either 
be empty or occupied by at most one calcium ion, the sublattice B take the values 0, 1 or 2 
corresponding to the closed, activating(open) and inhibiting state respectively. 

The dynamics of calcium ions in the model exhibit three stages. In the first, sponta- 
neous annihilation. If the site of the A sublattice is occupied then it becomes empty with 
probability g = (1 — p)a , here, p, related to the diffusion probability, a, related to the 
annihilation process. In the second, diffusion. One of the four nearest neighbor of site of 
the A sublattice, say site of sublattice B, is chosen at random. A calcium ion then hops 
from a site of one sublattice to a site of the other sublattice with probability p. In the third, 
Catalytic creation. One of the four nearest neighbor of site of the A sublattice, say site j of 
sublattice B, is chosen at random. If calcium channel j is open then a calcium ion is created 
at site i with probability r = (1 — — a). 

Coarse- Graining process — In the paper, neighboring microscopic sites q x q are grouped 
together into a CG cell, one can obtain a CG-lattice model with coarse cells. Fig. [T] shows 
an example of a coarse-graining lattice model with q = 2 and g = 4, denoted by solid square 
and dotted square respectively, here, q is even because of two kind of sublattice A and B. 
We define CG variables 

^ = E.,^, ^1 = E,,^, '^(^^- - ^2 = - 2) (1) 

Here, microscopic variable Si and xj denote the number of calcium ions at the A and B 
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FIG. 1: (Color online) Schematic illustration of coarse-graining lattice model. Circle sites denote 
sublattice A and square cites denote sublattice B. 

respectively, and satisfies the constraint < rj < < ai < < (72 < 

since each coarse cell contains microcells. Equivalently we may also consider the averaged 
version (termed below as coverage), s = {s = r], ai, a2). The dynamics of calcium ions 

on the CG model has also three processes: a. Spontaneous annihilation, b. Diffusion, c. 
Catalytic creation. The table [T] gives summary of processes and transition probability rates 
for CG-kMC. 



TABLE I: The processes and transition probability rates for CG-kMC. 



Process Change of coarse variables Coarse transition probability rate 
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III. RESULTS AND DISCUSSION 



Given a microscopic initial condition at random, following the aforementioned rules in 
table [U the coarse-grained calcium dynamics is computed with periodic boundary conditions. 
But it needs to make the computational demand of CG-kMC simulations per event the same 
as that of microscopic MC ones. 

We perform CG-kMC simulations and microscopic simulations on a square lattice with 
N X N = 200 X 200 sites, and plot the coverage 77, ai, a2 as a function of the parameter a 
in Fig. |2l where, ai denotes the density of open channels on sublattice B, 0-2 denotes the 
density of inhibited channels on sublattice B, and 77 denotes calcium ions on sublattice A. 
Firstly we notice that the coverage predicted from the CG-kMC simulations is in reasonably 
agreement with that of the microscopic MC ones. Excitedly, the CG-kMC predicts the phase 
transition point is in good agreement with that of microscopic MC simulations. Indeed, 
small quantitative differences near the critical point also exist, probably due in part to the 
fluctuation, but still relatively small. These findings validate the CG approach works well 
in simulating calcium dynamics. 




FIG. 2: CG-kMC simulation results for the coverage rj, cxi, C72 as a function of the parameter a for 
the case p = 0.5 and N x N = 200 x 200. The coverage vanish at the critical point Oc — 0.47. Lines 
indicate micro-simulations, symbols corresponds to CG ones with q = 40. 

To detect the phase transition point accurately, we need to choose appropriate size q x q 
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FIG. 3: (Color online) The deviation Aac of the phase transition points between CG-kMC and 
microscopic simulations vs the coarse proportion m for different system size N x N. 



of coarse cell. But what a suitable g is ? We define a coarse proportion m, m = N/q, namely, 
the square root of the number of coarse cells. Changing m we plot the deviation Aa^ of the 
phase transition points between CG-kMC and microscopic simulations for different system 
size N X N, a.s shown in Fig. [31 It can be seen that Aoc begin to decrease and then increase 
with the increment of m, especially, the minimal Aa^ occurs near the same point m = 6 for 
different N x N. When N x N = 720 x 720, there are two values of m corresponding to the 
minimal deviation, seeming to a small plain appears. Furthermore, the larger system size, 
the less deviation is. It is obvious that, there exists an appropriate range of m for coarse 
graining the system precisely. Therefore we can fix m and investigate the effects of system 
size on the phase transition point. 

For fixed m = 6, we plot the critical point ac as the function of system size N x N 
in Fig. m Apparently, ac increases monotonously as N x N increases, and approaches 
to 0.5. Theoretically speaking, this asymptotic value corresponds to the critical point of 
mean field (MF). The inset gives the scaling relation of the deviations and the system size, 
the scahng exponent is —0.531. To elucidate its accuracy, we have also carried out coarse 
grained simulations with m = 8 (not shown here), and obtained a similar asymptote and 
power law. Therefore, we can analyze the effects of the system size on the phase transition 
point according to the scaling law and detect the critical point accurately and rapidly. 

Finally, we exhibit the significant computational savings resulting from coarse-graining. 
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FIG. 4: The dependence of the critical point ac on the system size NxN for fixed m x m = 6 x 6. 
The inset gives the scahng law of the deviations Gcmf — clcCG and NxN, where Ocmf denotes the 
critical point of the MF model and OcCG denotes the critical point of the CG model, the scaling 
exponent is —0.531. 

as shown in Fig. [51 It can be seen that CPU time decreases monotonously as the size 
q X q of coarse cell increases. The larger size of coarse cell, the fewer CPU times. In this 
way, we can choose a bigger q to save computational time. In fact, the computer time in 
kMC simulation with global update, i.e., searching the entire lattice to identify the chosen 
site, scales approximately as 0(A^^), but O(m^) in CG-kMC simulation. Accordingly, a 
g-fold reduction in the number of sites results in reduced computer time by a factor of 
Therefore, coarse-graining can render MC simulation for the large scales feasible. 

IV. CONCLUSION 

In this paper, we proposed an extensive CG model that can properly describe calcium 
dynamics on ER membrane. By a great deal of computer simulations, we demonstrated 
our model is highly effective because the results of CG-kMC simulations are in very good 
agreement with that of MC ones for a wide range of model parameters. Interestingly, it was 
shown that there exists an appropriate range of coarse proportion m, corresponding to the 
best estimation on the phase transition point compared to the microscopic counterpart, and 
such m is almost insensitive to the change of the system size. This make it possible to select 
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FIG. 5: CPU time for CG-kMC as a function of coarse cell size qxq. 

an m without beforehand unwanted simulations for any real-world system size. Moreover, 
The CG-kMC method provides significant reduction in CPU while retaining very good ac- 
curacy in estimating the phase transition point. The larger the level of coarse-graining qxq 
is, the larger computational savings are, therefore we can obtain the phase transition point 
quickly. Using the CG model, we also found that the critical point increases monotonously 
as the system size increases. Especially, there exists a scaling relation between the deviations 
of the phase transition point and the system size. A major advantage of the coarse model 
is that they have a direct connection to the microscopic dynamics and can provide valuable 
insights. Due to its reasonable accuracy and low computational requirements, we anticipate 
that the methods outlined in this work for simple systems will find widespread use in many 
realistic systems. 

Acknowledgments 

This work was supported by the National Natural Science Foundation of China (Grant 
No. 11205002). C.S.S. was also supported by the Key Scientific Research Fund of Anhui 
Provincial Education Department (Grant No.KJ2012Al89). 




[1] M. J. Berridge, M. D. Bootman, and P. Lipp, Nature 395, 645 (1998). 



8 



[2] A. P. Thomas, G. S. J. Bird, G. Hajnoczky, L. D. Robb-Gaspers, and J. W. Putney, Faseb J. 
10, 1505 (1996). 

[3] J. J. Abramson, A. C. Zable, T. G. Favero, and G. Salama, J. Biol. Chem. 270, 29644 (1995). 
[4] J. H. Jaggar, V. A. Porter, W. J. Lederer, and M. T. Nelson, Am. J. Physiol. Cell Physiol. 

278, C235 (2000). 
[5] M. D. Stern, Biophys. J. 63, 497 (1992). 

[6] D. M. Bers, Excitation-Contraction Coupling and Cardiac Contractile Force, 2rd edn., 

Boston: Kluwer Academic Publishers, 215 (2001). 
[7] R. L. Winslow, A. Tanskanen, M. Chen, and J. L. Greenstein, Ann. N.Y. Acad. Sci. 1080, 

362 (2006). 

[8] R. Hinch, J. Greenstein, and R. Winslow, Prog. Biophys. Mol. Biol. 90, 136 (2006). 
[9] E. R. Higgins, P. Goel, J. L. Puglisi, D. M. Bers, M. Cannell, and J. Sneyd, J. Theor. Biol. 
247, 623 (2007). 

[10] R. Hinch, J. L. Greenstein, A. J. Tanskanen, L. Xu, and R. Winslow, Biophys. J. 87, 3723 
(2004). 

[11] M. Bar, M. Falcke, H. Levine, and L. S. Tsimring, Phys. Rev. Lett. 84, 5664 (2000). 
[12] J. W. Shuai and P. Jung, Proc. Natl. Acad. Sci. USA 100, 506 (2003). 

[13] M. A. Katsoulakis and D. G. Vlachos, J. Chcm. Phys. 119, 9412 (2003). 
[14] M. A. Katsoulakis, A. J. Majda, and D. G. Vlachos, Proc. Natl. Acad. Sci. USA 100, 782 
(2003). 

[15] M. A. Katsoulakis, A. J. Majda, and D. G. Vlachos, J. Comp. Phys. 186, 250 (2003). 
[16] A. Chatterjee and D. G. Vlachos, J. Comput.-Aided Mater. Des. 14, 253 (2007). 

[17] N. Guisoni and M. J. de Oliveira, Phys. Rev. E 71, 061910 (2005). 
[18] N. Guisoni and M. J. de Oliveira, Phys. Rev. E 74, 061905 (2006). 



9 



